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Sedimentation of a dispersed solid phase is widely encountered in applications and envi¬ 
ronmental flows, yet little is known about the behavior of finite-size particles in homo¬ 
geneous isotropic turbulence. To fill this gap, we perform Direct Numerical Simulations 
of sedimentation in quiescent and turbulent environments using an Immersed Boundary 
Method to account for the dispersed rigid spherical particles. The solid volume fractions 
considered are </) = 0.5 — 1%, while the solid to fluid density ratio Pp/pf = 1.02. The par¬ 
ticle radius is chosen to be approximately 6 Komlogorov lengthscales. The results show 
that the mean settling velocity is lower in an already turbulent flow than in a quiescent 
fluid. The reduction with respect to a single particle in quiescent fluid is about 12% and 
14% for the two volume fractions investigated. The probability density function of the 
particle velocity is almost Gaussian in a turbulent flow, whereas it displays large positive 
tails in quiescent fluid. These tails are associated to the intermittent fast sedimentation 
of particle pairs in drafting-kissing-tumbling motions. The particle lateral dispersion is 
higher in a turbulent flow, whereas the vertical one is, surprisingly, of comparable mag¬ 
nitude as a consequence of the highly intermittent behavior observed in the quiescent 
fluid. Using the concept of mean relative velocity we estimate the mean drag coefficient 
from empirical formulas and show that non stationary effects, related to vortex shedding, 
explain the increased reduction in mean settling velocity in a turbulent environment. 
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1. Introduction 

The gravity-driven motion of solid particles in a viscous fluid is a relevant process in a 
wide number of scientific and engineering applications ( Guazzelli fc Morris|2012 ). Among 
these we recall fluvial geomorphology, chemical engineering systems, as well as pollutant 
transport in underground water and settling of micro-organisms such as plankton. 

The general problem of sedimentation is very complex due to the high number of factors 
from which it depends. Sedimentation involves large numbers of particles settling in dif¬ 
ferent environments. The fluid in which the particles are suspended may be quiescent or 
turbulent. Particles may differ in size, shape, density and stiffness. The range of spatial 
and temporal scales involved is wide and the global properties of these suspensions can 
be substantially different from one case to another. Because of these complexities, our 
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general understanding of the problem is still incomplete. 


1.1. Settling in a quiescent fluid 

One of the earliest investigations on the subject at hand is Stokes’ analysis of the sedi¬ 
mentation of a single rigid sphere through an unbounded quiescent viscous fluid at zero 
Reynolds number. This led to the well-known formula that links the settling velocity 
to the sphere radius, the solid to fluid density ratio and the viscosity of the fluid that 
bears his name. Later, the problem was studied both theoretically and experimentally. 


Hasimoto (1959) obtained expressions for the drag force exerted by the fluid on three 


different cubic arrays of rigid spheres. These relate the drag force only to the solid volume 
fraction, but were derived under the assumption of very dilute suspensions and Stokes 


flow. The formulae have later been revisited by 

Sangani & Acrivos 

(1982 

). A different 

approach was instead pursued by 

Batchelor 

(1972 

) who found a relation between the 


mean settling velocity and the solid volume fraction by using conditional probability 
arguments. When the Reynolds number of the settling particles {Ret) becomes finite, 
the assumption of Stokes flow is less acceptable (especially for Ret > 1). The fore-aft 
symmetry of the fluid flow around the particles is broken and wakes form behind them. 
Solutions should be derived using the Navier-Stokes equations, but the nonlinearity of 
the inertial term makes the analytical treatment of such problems extremely difficult. 
For this reason theoretical investigations have progressively given way to experimental 
and numerical approaches. 

The first remarkable experimental results obtained for creeping flow were those by [Richard - 
son & Zaki (19541. These authors proposed an empirical formula relating the mean set¬ 


tling velocity of a suspension to its volume fraction and to the settling velocity of an 
isolated particle. This formula is believed to be accurate also for concentrated suspen¬ 
sions (up to a volume fraction (j) of about 25%) and for low Reynolds numbers. Subsequent 
investigations improved the formula so that it could also be applied in the intermediate 
Reynolds numbers regime ( Garside fc Al-Dibouni|[l977 Di Felice||1999 1. 

Efficient algorithms and sufficient computational power have become available only rel¬ 
atively recently and since then many different numerical methods have been used to 


improve our understanding of the problem (Prosperetti 2015). Among others we recall 


the dynamical simulations performed by Ladd (1993), the finite-elements simulations 


of Johnson & Tezduyar (1996), the force-coupling method simulations by Climent & 


Maxey (2003), the Lattice-Boltzmann simulations of Yin & Koch (2007), the Oseenlet 
simulations by Pignatel et al. (20111, and the Immersed Boundary simulations of Kempe 


& Frohlich (2012) and Uhlmann & Doychev (2014). Thanks to the most recent tech¬ 


niques it has become feasible to gain more insight on the interactions among the different 


phases and the resulting microstructure of the sedimenting suspension (Yin & Koch 2007 


Uhlmann & Doychev 2014). Uhlmann & Doychev (2014), most recently, simulated the 


settling of dilute suspensions with particle Reynolds numbers in the range 140 — 260 and 
studied the effect of the Archimedes number (namely the ratio between gravitational and 
viscous forces) on the microscopic and macroscopic properties of the suspension. These 
authors observe an increase of the settling velocity at higher Archimedes, owing to parti¬ 
cle clustering in a regime where the flow undergoes a steady bifurcation to an asymmetric 
wake. Settling in stratified environments has also been investigated experimentally, i.e. 


by Bush et al. (2003), and numerically, i.e. by Doostmohammadi & Ardekani (2015). 
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The investigations previously reported consider the settling of particles in quiescent or 
uniform flows. There are many situations though, where the ambient fluid is in fact 
nonuniform or turbulent. As in the previous case, the first approach to this problem was 


analytical. In the late 40’s and 50’s Tchen (1947) and later Corrsin & Lumley (1956) 


proposed an equation for the motion of a small rigid sphere settling in a nonuniform 
flow. In the derivation, they assumed the particle Reynolds number to be very low so 
that the viscous Stokes drag for a sphere could be applied. The added mass and the 
augmented viscous drag due to a Basset history term were also included. |Maxey fc Riley] 
(1983) corrected these equations including also the Faxen forces due to the unsteady 


Stokes flow. 

In a turbulent flow many different spatial and temporal scales are active. Therefore the 
behaviour and motion of one single particle does not depend only on its dimensions and 
characteristic response time, but also on the ratios among these and the characteristic 
turbulent length and time scales. The turbulent quantities usually considered are the Kol¬ 
mogorov length and time scales which are related to the smallest eddies. Alternatively, 
the integral lengthscale and the eddy turnover time can also be used. It is clear that a 
particle smaller than the Kolmogorov legthscale will behave differently than a particle of 
size comparable to the energetic flow structures. A sufficiently large particle with a char¬ 
acteristic time scale larger than the timescale of the velocity fluctuations will definitely 
be affected by and affect the turbulence. A smaller particle with a shorter relaxation time 
will more closely follow the turbulent fluctuations. When particle suspensions are con¬ 
sidered, the situation becomes even more complicated. If the particles are solid, smaller 
than the Kolmogorov lengthscale and dilute, the turbulent flow field is unaltered (i.e., 
one-way coupling). Interestingly, the turbulent dynamics is instead altered by microbub¬ 
bles. The presence of these microbubbles leads to relevant drag reduction in boundary 
layers and shears flows (e.g. Taylor-Couette flow) ( Sugiyama et aZ.||2008| Ceccio|[2M0 ). If 
the mass of the dispersed phase is similar to that of the carrier phase, the influence of the 
solid phase on the fluid phase cannot be ignored (i.e., two-way coupling). Interactions 
among particles (such as collisions) must also be considered in concentrated suspensions. 
This last regime is described as four-way coupling (|Elgobashi||1991| |Balachandar & Eaton 


Because of the difficulty of treating the problem analytically, the investigations of 
the last three decades have mostly been either experimental or numerical. In most of 
the numerical studies heavy and small particles were considered. The reader is referred 


to 

Toschi & Bodenschatz 

(2009 

) for 

a more detailed review than the short summary 

reported here. 

Wang &: Maxey 

(1993 

) studied the settling of dilute heavy particles in 


homogeneous isotropic turbulence. The particle Reynolds number based on the relative 
velocity was assumed to be much less than unity so that Stokes drag force could be used 
to determine the particle motion. These authors show that heavy particles smaller than 
the Kolmogorov lengthscale tend to move outward from the center of eddies and are often 
swept into regions of downdrafts (the so called preferential sweeping later renamed fast- 
tracking) . In doing so, the particle mean settling velocity is increased with respect to that 
in a quiescent fluid. A series of studies confirmed and extended these results examining 


(Aliseda et al. 


(Siewert et al. 


particle clustering (Bee et al. 2014 Gustavsson et al. 2014), preferential concentration 


2002), the effects of the particle shape, orientation and collision rates 


2014), as well the effects of one- or two-way coupling algorithms (Bosse 


et a/.|2006 1, to mention few aspects. Numerous experimental studies were also performed 
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in order to confirm these results and to study the turbulence modulation due to the 
presence of particles ( Hwang fc Eaton|[2006 1. 

The results on the mean settling velocities of particles of the order or larger than the 
Kolmogorov scale are not conclusive. Good et al.\ ( 2014| studied particles smaller than 
the Kolmogorov scale and with density ratio 0(1000), whereas Variano (experiments; 


private communication) and Byron (20151 studied finite-size particles at density ratios 


comparable to ours. Good et al. (2014) found that the mean settling velocity is reduced 


only when nonlinear drag corrections are considered in a one-way coupling approach 
when particles have a long relaxation time (a linear drag force would always lead to a 
settling velocity enhancement). For finite-size almost neutrally-buoyant particles, Vari¬ 
ano and Byron (2015) observe instead that the mean settling velocity is smaller than in 


a quiescent fluid. In relative terms, the settling velocity decreases more and more as the 
ratio between the turbulence fluctuations and the terminal velocity of a single particle 
in a quiescent fluid increases. It is generally believed that the reduction of settling speed 
is due to the non-linear relation between the particle drag and the Reynolds number. 
Nonetheless, unsteady and history effects may also play a key role (Olivieri et aL|[2M4 


Bergougnoux et al. 2014). Tunstall & Houghton (1968) demonstrated already in 1968 


that the average settling velocity is reduced in a flow oscillating about a zero mean, due 


to the interactions of the particle inertia with a non-linear drag force. Stout et al. (1995) 


tried to motivate these findings in terms of the relative motion between the fluid and the 
particles. When the period of the fluid velocity fluctuations is smaller than the particle 
response time, a significant relative motion is generated between the two phases. Due 
to the drag non-linearity, appreciable upward forces can be produced on the particles 
thereby reducing the mean settling velocity. 

Unsteady effects may become important when considering suspensions with moderate 


particle-fluid density ratios, as suggested by Mordant & Pinton (2000) and Sobral, 


Oliveira & Cunha (20071. The former studied experimentally the motion of a solid sphere 


settling in a quiescent fluid and explain the transitory oscillations of the settling velocity 
found at Re « 0(100) by the presence of a transient vortex shedding in the particle 
wake. The latter, instead, analyzed an equation similar to that proposed by [Maxey ^ 
Riley (1983), and suggested that unsteady hydrodynamic drags might become important 


when the density ratio approaches unity. 

1.3. Fully resolved simulations 

As already mentioned, most of the numerical studies of settling in turbulent flows used ei¬ 
ther one or two-way coupling algorithms. In order to properly understand the microscop¬ 
ical phenomena at play, it would be ideal to use fully resolved simulations. An algorithm 
often used to accomplish this is the Immersed Boundary Method with direct forcing for 


the simulation of particulate flows originally developed by Uhlmann (2005). The code 


was later used to study the clustering of non-colloidal particles settling in a quiescent 


environment (Uhlmann & Doychev 2014). With a similar method Lucci et al. (2010) 


studied the modulation of isotropic turbulence by particles of Taylor length-scale size. 


Recently, Homann et al. (2013) used an Immersed Boundary Fourier-spectral method 


to study finite-size effects on the dynamics of single particles in turbulent flows. These 
authors found that the drag force on a particle suspended in a turbulent flow increases as 
a function of the turbulent intensity and the particle Reynolds number. We recently used 


a similar algorithm to examine turbulent channel flows of particle suspensions (Picano, 
Breugem fc Brandt||2015 ). 


The aim of the present study is to simulate the sedimentation of a suspension of parti¬ 
cles larger than the Kolmogorov lengthscale in homogeneous isotropic turbulence with a 
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finite difference Immersed Boundary Method. We focus on particles slightly denser than 
the suspending fluid {pp/pf = 1.02) and investigate particle and fluid velocity statistics, 
non-linear and unsteady contributions to the overall drag and turbulence modulation. 
The suspensions considered in this work are dilute = 0.5—1%) and monodispersed. The 
same simulations are also performed in absence of turbulence to appreciate differences of 
the particle velocity statistics in the two different environments. Due to the size of the 
particles considered it has been necessary to consider very long computational domains 
in the settling direction, especially for the quiescent environment. In the turbulent cases, 
smaller domains provide converged statistics since the particle wakes are disrupted faster. 
The parameters of the simulations have been inspired by the experiments by Variano, 


Byron (2015) and co-workers at UC Berkeley. These authors investigate Taylor-scale 


particles in turbulent aquatic environments using Refractive-Index-Matched Hydrogel 
particles to measure particle linear and angular velocities. 

Our results show that the mean settling velocity is lower in an already turbulent flow 
than in a quiescent fluid. The reduction is about 12% and 14% for the two volume 
fractions investigated. By looking at probability density functions pdf of the settling 
velocities, we observe that the pdf is well approximated by a Gaussian function centered 
around the mean in the turbulent cases. In the laminar case instead, the pdf shows 
a smaller variance and a larger skewness, indicating that it is more probable to find 
particles settling more rapidly than the mean value rather than more slowly. These 
events are associated to particle-particle interactions, in particular to the drifting-kissing- 
tumbling motion of particle pairs. We also calculate mean relative velocity fields and 
notice that vortex shedding occurs around each particle in a turbulent environment. 
Using the concept of mean relative velocity we calculate a local Reynolds number and 
the mean drag coefficient from empirical formulas to quantify the importance of unsteady 
and history effects on the overall drag, thereby explaining the reduction in mean settling 
velocity. In fact, these terms become important only in a turbulent environment. 


2. Methodology 

2.1. Numerical Algorithm 

Different methods have been proposed in the last years to perform Direct Numerical Sim¬ 
ulations of multiphase flows. The Lagrangian-Eulerian algorithms are believed to be the 
most appropriate for solid-fluid suspensions (Ladd fc Verberg|2001 Zhang & Prosperetti 


[2010 Lucci et al. 2010 Uhlmann & Doychev 2014). In the present study, simulations have 
been performed using a tri-periodic version of the numerical code originally developed 


by Breugem (2012) that models the coupling between the solid and fluid phases. The 


Eulerian fluid phase is evolved according to the incompressible Navier-Stokes equations, 

V-u/ = 0 (2.1) 


dw 1 

—^ -I- Uf • Vuf =-Vp -I- vW^Vif -b f (2.2) 

at ^ Pf 

where Uf,pf and v = p/pf are the fluid velocity, density and kinematic viscosity respec¬ 
tively {p is the dynamic viscosity), while p and f are the pressure and the force field used 
to mantain turbulence and model the presence of particles. The particles centroid linear 
and angular velocities, Up and ujp are instead governed by the Newton-Euler Lagrangian 
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equations, 


= Pf T • n dS' + (pp - pf) Vpg 


IdVj 

_ f r X T ■ ndS 
^ dt Jav^ 


(2.3) 

(2.4) 


where Vp = 47Ta^/3 and Ip = 2ppVpC?jh are the particle volume and moment of inertia, 
with a the particle radius; g is the gravitational acceleration; r = —pi + 2/iE is the fluid 
stress, with I the identity matrix and E = /2 the deformation tensor; r 

is the distance vector from the center of the sphere while n is the unit vector normal to 
the particle surface dVp. Dirichlet boundary conditions for the fluid phase are enforced 
on the particle surfaces as u/|aVp — Up + ujp x r. 

In the numerical code the coupling between the solid and fluid phases is obtained by 
using an Immersed Boundary Method. The boundary condition at the moving particle 
surface (i.e. u/|aVp = Up + Wp x r) is modeled by adding a force field on the right-hand 
side of the Navier-Stokes equations. The problem of re-meshing is therefore avoided and 
the fluid phase is evolved in the whole computational domain using a second order finite 
difference scheme on a staggered mesh. The time integration is performed by a third order 
Runge-Kutta scheme combined with a pressure-correction method at each sub-step. The 


same integration scheme is also used for the Lagrangian evolution of eqs. (2.3) and (2.4). 


The forces exchanged by the fluid and the particles are imposed on Nl Lagrangian 
points uniformly distributed on the particle surface. The force F; acting on the I — 
th Lagrangian point is related to the Eulerian force field f by the expression f(x) = 
SzSl F;(5d(x — Xi)AV;. In the latter AV; represents the volume of the cell containing the 
I — th Lagrangian point while Sd is the Dirac delta. This force field is obtained through 
an iterative algorithm that mantains second order global accuracy in space. Using this 


IBM force field eqs. (2.3) and (2.4) are rearranged as follows to maintain accuracy, 


7 ^ r! r 

Pp'^p^ = Pf^ ^fdV+{pp-pf)Vpg (2.5) 

( 2 . 6 ) 


/ = 1 
Ni 


dUJp T-tATX ^ f 

= + r^^,dV 

1=1 U 


where the second terms on the right-hand sides are corrections to account for the inertia 


of the fictitious fluid contained within the particle volume. In eqs. (2.5),(2.6) r; is simply 


the distance from the center of a particle. Particle-particle interactions are also consid¬ 
ered. When the gap distance between two particles is smaller than twice the mesh size, 


lubrication models based on Brenner’s asymptotic solution (Brenner 1961) are used to 


correctly reproduce the interaction between the particles. A soft-collision model is used 
to account for collisions among particles with an almost elastic rebound (the restitution 
coefficient is 0.97). These lubrication and collision forces are added to the right-hand side 
of eq. (2.5). More details and validations of the numerical code can be found in |Breug^ 
(2012), Lambert et al. (2013), Lashgari et al. (2014) and Picano et al. (2015). 


2.2. Parameter setting 

Sedimentation of dilute particle suspensions is considered in an unbounded computational 
domain with periodic boundary conditions in the x, y and z directions. Gravity is chosen 
to act in the positive z direction. A zero mass flux is imposed in the simulations. A 
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cubic mesh is used to discretise the computational domain with eight points per particle 
radius, a. Non-Brownian rigid spheres are considered with solid to fluid density ratio 
Pp!Pf ~ 1-02. Hence, we consider particles slightly heavier than the suspending fluid. Two 
different solid volume fractions (j) = 0.5% and 1% are considered. In addition to the solid 
to fluid density ratio Pp/pf and the solid volume fraction (f, it is necessary to introduce 
another nondimensional number. This is the Archimedes number (or alternatively the 
Galileo number Ga = -s/Ar), 


Ar = 


{jj - 1) g(2«)^ 


(2.7) 


a nondimensional number that quantifies the importance of the gravitational forces acting 
on the particle with respect to viscous forces. In the present case the Archimedes number 
Ar = 21000. Using the particle terminal velocity Vt we define the Reynolds number 


Ret = 2avt/v. This can be related by empirical relations to the drag coefficient of an 
isolated sphere when varying the Archimedes number, Ar. Different versions of these 
empirical relations giving the drag coefficient as a function of Ret and Ar have been 

proposed. As Yin & Koch 

(2007) we will use the following relations. 

( 24 

= i S: ■ 

1 Ret . 

since C/j = 4Ar/(3Ret) ( 

Ar=l 

[ 18i?et 

1 + 0.1315Ref ^ Q ^ 20 ^ 

L + 0.1935i?e?'6305] ^ 20 < Ret < 260 

Yin & Koch||2007 ) we finally write, 

1 + 0.1315i?e^° ®^”° °3'°®'“^'=‘^ , 0.01 < Ret ^ 20 g. 

1 + 0.1935i?e° ®305] ^ 20 < Ret < 260 


The Reynolds number calculated from eq. (2.9) is approximately 188 for Ar = 21000. 


In order to generate and sustain an isotropic and homogeneous turbulent flow field, a 
random forcing is applied to the first s hell of wave vectors. We consider a (5-correlate d 


in time forcing of fixed amplitude /o (Vincent & Meneguzzi 


1991 


Zhan et al. 


2014). 


The turbulent field, alone, is characterized by a Reynolds number based on the Taylor 
microscale, Re\ = Xu'/v, where u' is the fluctuating velocity and A = \/Vovu'"^ je is the 
transverse Taylor length scale. This is about 90 in our simulations. The ratio between 
the grid spacing and the Kolmogorov lengthscale ry = (y^ jeY^'^ (where e is the energy 
dissipation) is approximately 1.3 while the particle diameter is circa I 277 . Finally, the 
ratio among the expected mean settling velocity and the turbulent velocity fluctuations 
is Vtlu! = 3.4. The parameters of the turbulent flow field are summarized in table For 
the definition of these parameters, the reader is referred to Pope (2000). 


2.3. Validation 

To check the validity of our approach we performed simulations of a single sphere settling 
in a cubic lattice in boxes of different sizes. Since this is equivalent to changing the solid 
volume fraction, we compared our results to the analytical formula derived by |Hasimoto 
(1959) and Sangani & Acrivos ( |1982 ), 


|Vt| = ^ = 1 - 1.7601.^1/3 


( 2 . 10 ) 


where 


|V,| = 






( 2 . 11 ) 
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Figure 1. Terminal velocity of a periodic array of spheres. Present results Ret = 1 denoted by 
red circles are compared with the ones obtained by Yin and K och by black squares at the same 
Ret and with the analytical solution for Ret = 0 of eq. (2.101. 


rj/{2a) u' k \/{2a) Rex e Te ReLg 

0.084 0.30 0.13 1.56 90 0.0028 46.86 1205 

Table 1. Turbulent flow parameters in particle units, where k is the turbulent kinetic energy, A 
is the Taylor microscale, Te = k/e is the eddy turnover time and Rerg is the Reynolds number 
based on the integral lengthscale Lq = 


is t he St okes settling velocity. In terms of the size of the computational domain, 1^, 
eq. (2.10) can also be rewritten as V* = 1 — 1.7601(2a/Zz)(7T/6)^/^. In figure in we show 
the results obtained with our code for Ret = 1 together with the data by [Yin fc Koch| 
(2007) and the analytical solution. Although the analytical solution was derived with the 
assumption of vanishing Reynolds numbers, we find good agreement among the various 
results. 

The actual problem arises when considering particles settling at relatively high Reynolds 
numbers. If the computational box is not sufficiently long in the gravity direction, a 
particle would fall inside its own wake (due to periodic boundary conditions), thereby 
accelerating unrealistically. Various simulations of a single particle falling in boxes of dif¬ 
ferent size were preliminarily carried out, in particular 48a x 48a x 48a, 32a x 32o x 96a 
and 32a x 32a x 320a. The first two boxes turn out to be unsuitable for our purposes. 
We find a terminal Reynolds number Ret = 200 in the longest domain considered, which 
corresponds to a difference of about 6% with respect to the value obtained from the 
empirical relations (2.9). As reference velocity we use the value obtained from simula¬ 
tions performed in the largest box at a solid volume fraction tw o orders of magnitude 
smaller than the cases under investigation, (j> = 5 ■ 10“® (as in Uhlmann & Doychev 


2014), corresponding to a terminal velocity such that Ret = 195, 4% larger than the 


value from the empirical relations (2.9). Further increasing the length in the 2 ; direction 
would make the simulations prohibitive. Note also that simulations in a turbulent envi¬ 
ronment turn out to be less demanding as turbulence disrupts and decorrelates the flow 
structures induced by the particles. The final choice was therefore a computational box 
of size 32o x 32a x 320a with 256 x 256 x 2560 grid points, 391 particles for (j) = 0.5% 
and 782 particles for (p = 1%. In all cases, the particles are initially distributed randomly 
in the computational volume with zero initial velocity and rotation. 
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Figure 2. Instantaneous snapshot of the velocity component perpendicular to gravity on dif¬ 
ferent orthogonal planes, together with the corresponding particle position for f = 0.005. A 
zoomed view of a particular section is also shown. 

A snapshot of the suspension flow for cf = 0.5% is shown in figureThe instantaneous 
velocity component perpendicular to gravity is shown on different orthogonal planes. 

The simulations were run on a Cray XE6 system at the PDC Center for High Perfor¬ 
mance Computing at the KTH, Royal Institue of Technology. The fluid phase is evolved 
for approximately 6 eddy turnover times before adding the solid phase. The simulations 
for each solid volume fraction were performed for both quiescent fluid and turbulent flow 
cases in order to compare the results. Statistics are collected after an initial transient 
phase of about 4 eddy turnover times for the turbulent case and 15 relaxation times 
(Tp = 2ppaf/{9pfv)) for the quiescent case. Defining as reference time the time it takes 
for an isolated particle to fall over a distance equal to its diameter, 2ajvt^ the initial 
transient corresponds to approximately 170 units. Statistics are collected over a time 
interval of 500 and 300 in units of 2a/vt for the quiescent and turbulent cases, respec¬ 
tively. Differences between the statistics presented here and those computed from half 
the samples is below 1% for the first and second moments. 


3. Results 

We investigate and compare the behavior of a suspension of buoyant particles in a 
quiescent and turbulent environment. The behaviour of a suspension in a turbulent flow 
depends on both the particles and turbulence characteristic time- and length-scales: ho¬ 
mogeneous and isotropic turbulence is defined by the dissipative, Taylor and integral 
scales, whereas the particles are characterized by their settling velocity Vt and by their 
Stokesian relaxation time Tp = 2ppof /{9pfi') = 11.1 (time is scaled by {2a/vt) through¬ 
out the paper). A comparison between characteristic time scales is given by the Stokes 
number, i.e. the ratio between the particle relaxation time and a typical flow time scale 
Stf = Tp/Tf. In the present cases, the Stokes number based on the dissipative scales 
(time and velocity) is Stp = Tp/Tx = 8.1 so the particles are inertial on this scale. In 
addition, because the particles are about 12 times larger than the Kolmogorov length 
and fall about 16 times faster than the Kolmogorov velocity scale, we can expect that 
motions at the smallest scales weakly affect the particle dynamics. 

Considering therefore the large-scale motions, we introduce an integral-scale Stokes 
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Figure 3. Probability density functions of the particles velocities along the directions a) parallel 
Vp^T and b) perpendicular Vp^n to gravity for (j) = 0.5% and 1%. The velocities are normalized 
by vt- Quiescent fluid cases denoted by blue curve with squares for 0 = 0.5% and green curve 
with triangles for (j) = 1%, turbulent flow data by red curve with circles for (p = 0.5% and black 
curve with downward-pointing triangles for 0 = 1%. In the insets we show the pdfs in lin-log 
scale. 


number = Tp/Tg = 0.24. This value of reflects the fact that the particles are 
about 20 times smaller than the integral length scale Lq. The strong coupling between 
particle dynamics and turbulent flow held occurs at scales of the order of the Taylor 
scale for the present cases. Indeed, the Taylor Stokes number is St\ = Tp/T\ = 2.1 
with T\ = A/m'. It should be noted that the Taylor length is slightly larger than the 
particle size, 3.1o, while particles fall around 3.4 times faster than typical huid velocity 
fluctuations, Vt/u' = 3.4. Hence particles are strongly inhuenced by the fluid huctuations 
occurring at scales of the order of A. 


3.1. Particle statistics 

We start by comparing the single-point how and particle velocity statistics for the two 
cases studied, i.e. quiescent and turbulent how. The results are collected when a statisti¬ 
cally steady state is reached. Due to the axial symmetry with respect to the direction of 
gravity, we consider only two velocity components for both phases, the components par¬ 
allel and perpendicular to gravity, Va,T and Va,n respectively, where a = f, p indicates 
the solid and huid phases. 

In hgurej^we report the probability density function of the particle velocities for both 
volume fraction investigated here, (j) = 0.5% and 1%; the moments extracted from these 
distribution are summarized in table The data in [^) show the pdf for the component 
of the velocity aligned with gravity Vp^r normalized by the settling velocity for ^ 0 

in a quiescent environment, Vt- This is extracted from the simulation of a very dilute 
suspension discussed in section [T3| 

In the quiescent cases, the mean settling velocity slightly reduces when increasing the 


volume fraction cp, in agreement with the hndings of Richardson & Zaki (1954) and 


Di Felice (19991 among others. The sedimentation velocity decreases to 0.96 at </ = 0.5% 


and 0.93 at (/ = 1%. Di Felice (19991 investigated experimentally the settling velocity of 


dilute suspensions of spheres ((/ = 0.5%) with density ratio 1.2 in quiescent fluid, for a 
large range of terminal Reynolds numbers (from 0.01 to 1000). Following the empirical 


fit proposed in Di Felice (19991, we obtain (V/,t) — 0.98 at (/ = 0.5%, approximately 1% 
larger than our result. On the other hand, the formula suggested for the intermediate 
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regime in Di Felice (19991 and Yin & Koch (2007), Ret < 150, would give an estimated 


value of 0.88, 6% lower than our result. 

The mean settling velocity for the quiescent case at </> = 0.5% is instead close to the 


estimated value from eq. (2.9). This is in agreement with what reported in Uhlmann & 


Doychev (2014). These authors found that the particle mean settling velocity of a dilute 


suspension increases above the reference value only when the Archimedes number Ar is 
larger than approximately 24000 and clustering occurs. In our case Ar ~ 21000 and no 
clustering is noticed accordingly with their findings. 

Interestingly, we observe an additional non-negligible reduction when a turbulent back¬ 
ground flow is considered, in our opinion the main result of the paper. The reduction of 
the mean settling velocity (Ip,r) is 12% at ^ = 0.5% and 14% at ^ = 1%, see table 
This result unequivocally shows that the turbulence reduces the settling velocity of a 
suspension of finite-size buoyant particles, in agreement with the experimental findings 
by Variano and co-workers (Variano, private commnnication) and Byron (2015). We also 
note that the reduction of the settling velocity with the volume fraction is less important 
for the turbulent cases. 

Looking more carefully at the pdfs we note that fluctuations are, as expected, larger 
in a turbulent environment. In addition, the vertical particle velocity fluctuations are 
the largest component in a quiescent fluid, whereas in a turbulent flow the fluctuations 
are largest in the horizontal directions, as summarized in table In the quiescent case 
the rms of the tangential velocity ^ is 0.15 and 0.17 for (j) = 0.5% and (f = 1% 
respectively, while it increases up to 0.23 in the corresponding turbulent cases. The 
difference in the width of the pdf is particularly large in the directions normal to gravity 
where the rms of the variance cry^ ^ is 0.3 for both turbulent cases, while it is 5 and 
4 times smaller for the quiescent flows a.t (j) = 0.5% and (j) = 1%. We believe that the 
interactions among the particle wakes, mainly occurring in the settling direction, promote 
the higher vertical velocity fluctuations found in the quiescent cases. The shape of the pdf 
is essentially Gaussian for the turbulent cases, showing vanishing skewness and normal 
flatness. Interestingly an intermittent and skewed behaviour is exhibited in the quiescent 
cases. The flatness K is around 9 for both components at ^ = 0.5% and slightly reduces 
to 5.5 at (j) = 1%. The settling velocity of the quiescent cases is skewed towards intense 
fluctuations in the direction of the gravity. The skewness S is higher for the more dilute 
case, being 1.26 at ^ = 0.5% and 0.7 at ^ = 1%. 

We interpret the intermittent behavior suggested by values of AT > 3 by the collective 
dynamics of the particle suspension. The significant tails of the pdfs shown in figure [^) 
are indeed associated to a specific behavior: as particles fall, they tend to be accelerated 
by the wakes of other particles, before showing drafting-kissing-tumbling behavior (|Fortes| 


et al. 1987). Snapshots of the drafting-kissing-tumbling behavior between two spherical 


particles are shown in figure]^ When this close interaction occurs, particles are found to 
fall with velocities that can be more than twice the mean settling velocity In the 

quiescent case the fluid is still, the wakes are the only perturbation present in the field 
and are long and intense so their effect can be felt far away from the reference particle. 
The more dilute the suspension the more intermittent the particle velocities are. On the 
contrary, when the flow is turbulent the wakes are disrupted quickly and therefore fewer 
particles feel the presence of a wake. The velocity disturbances are now mainly due to 
turbulent eddies of different size that interact with the particles to increase the variance 
of the velocity homogeneously along all directions leading to the almost perfect normal 
distributions shown above, with variances similar to those of the turbulent fluctuations. 
The features of the particle wakes will be further discussed in this manuscript to support 
the present explanation. 
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{Vp.z) 

Quiescent f = 0.5% 
-tO.96 
-t0.15 
-tl.26 
-t9.65 

Turbulent f = 0.5% 
+0.88 
+0.23 
+0.01 
+2.92 

Quiescent 0=1% 
+0.93 
+0.17 
+0.70 
+6.01 

Turbulent 0=1% 
+0.86 
+0.23 
+0.01 
+3.15 

{Vp,n) 

-tO.33- lO""* 

-1.93 • 10"® 

-8.78- 10"^ 

-0.97- 10"® 

e'Vp,„ 

-t0.06 

+0.31 

+0.08 

+0.31 


-1.22 • 10"® 

+0.04 

-3.17- 10"® 

+0.07 


-t8.95 

+2.78 

+5.55 

+2.80 


Table 2. First four central moments of the probability density functions of Vp^r and Vp," 
normalized by vt- Svp,^ (>S'vp,„) and are respectively the skewness and the flatness 

of the pdfs. 


We also note that the sedimenting speed in the quiescent fluid is determined by two 
opposite contributions. The excluded volume effect that contributes to a reduction of the 
mean settling velocity with respect to an isolated sphere and the pairwise interactions 
(the drafting-kissing-tumbling), increasing the mean velocity of the two particles involved 
in the encounter. To try to identify the importance of the drafting-kissing-tumbling effect, 
we fit the left part (where no intermittent behaviour is found) of the pdf pertaining the 
quiescent case at 4> = 0.5% with a gaussian function. The mean of Vp^r is reduced to 
about 0.93 (value due only to the hindrance effect) instead of 0.96 in the full simulation; 
thereby the increment in mean settling velocity due to drafting-kissing-tumbling can be 
estimated to be of about 3%. 

The probability density function of the particle angular velocities are also different in 
quiescent and turbulent flows. These are shown in figures [^) and b) for rotations about 

an axis parallel to gravity, \ajp^p\ = and orthogonal to it, settling 

direction the peak of the pdf is always at \oJp^z\ = 0. As for the translational velocities, 
the pdfs are broader in the turbulent cases. Due to the interaction with turbulent eddies, 
particles tend also to spin faster around axis perpendicular to gravity. From figure [^) 
we see that the modal value slightly increases in the quiescent cases when increasing 
the volume fraction. In the turbulent cases, the modal value is more than 3 times the 
values of the quiescent cases and the variance is also increased. Unlike the quiescent 
cases, the curves almost perfectly overlap for the two different (f, meaning that turbulent 
fluctuations dominate the particle dynamics. Turbulence hinders particle hydrodynamic 
interactions. 

Figure shows the temporal correlations of the particle velocity fluctuations, 


\^'^) — 2 

^P,T 

Z3 ^ _ ApAP’ t + At)) 

^Vn^nX^^) — 2 

^ V 

^p,n 


(3.1) 

(3.2) 


for the turbulent and quiescent cases at f) = 0.5% and (j) = 1%. Focusing on the data 
at the lower volume fraction, we observe that the particle settling velocity decorrelates 
much faster in the turbulent environment, within At ^ 50, while it takes around one 
order of magnitude longer in a quiescent fluid. Falling particles may encounter intense 
vortical structures that change their settling velocity. The turbulence strongly alters 
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Figure 4. Drafting-kissing-tumbling behavior among two spherical particles in the quiescent 
case with f = 0.5%. The particles are coloured with the absolute value of their velocity compo¬ 
nent in the direction of gravity. Three particles are labeled with ”A,B,C” in order to show how 
accelerated the two interacting particles are compared to the others. 
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Figure 5. Probability density functions of a) and b) + ^p,y for cj) = 0.5% and 

1%. The angular velocities are normalized by vt/(2a), the settling velocity of a single particle 
in a quiescent environment and its diameter. Quiescent fluid cases denoted by blue curve with 
squares for (j) = 0.5% and green curve with triangles for (j) = 1%, turbulent flow data by red 
curve with circles for (j) = 0.5% and black curve with downward-pointing triangles for tj) = 1%- 






At Vj/(2a) At Vj/(2a) 

Figure 6. Time correlations of the particle velocity fluctuations. The blue curve with circles 
represents the correlation of (Rv^v^), while the red curve with diamonds is used for the 
correlation of Vp^„ {Bv„v„)- a) Quiescent case with (j) = 0.5%, b) Turbulent case with cji = 0.5%, 
c) Quiescent case with (j) = 1.0% and d) Turbulent case with (j) = 1.0%. 


the fluid velocity field seen by the particles, which in the quiescent environment is only 
constituted by coherent long particle wakes. This results in a faster decorrelation of the 
velocity fluctuations along the settling direction in the turbulent environment. Moreover, 
Rv„v„ crosses the null value earlier than for the settling velocity component. This result 
is not surprising since the particle wakes develop only in the settling direction. 
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The normal velocity correlation Rv„v„ of the turbulent case oscillates around zero 
before vanishing at longer times. We attribute this to the presence of the large-scale 
turbulent eddies. As a settling particle encounters sufficiently strong and large eddies, 
its trajectory is swept on planes normal to gravity in an oscillatory way. To provide an 
approximate estimate of this effect, we consider as a first approximation the turbulent 
flow seen by the particles as frozen since the particles fall at a higher velocity than 
the turbulent fluctuations (3.4 times). Since the strongest eddies are of the order of 
the transversal integral scale we can presume that these structures are responsible for 
this behavior. In particular, the transversal integral scale Lt = To/2 — 8 and u'^^ns = 


0.3, so we expect a typical period oi t = Lt/u'^ 
oscillations found for both turbulent cases, i.e. t ? 


.jg ~ 26 which is of the order of the 
20. Note that a similar behavior has 
been observed byjWang fc Maxey|(|1993[) for sufficiently small and heavy particles, termed 


the preferential sweeping phenomenon. 

The same process can be interpreted in terms of crossing trajectories and continuity 


effects as described by Csanady (1963). An inertial particle falling in a turbulent envi¬ 


ronment changes continuously its fluid-particle neighborhood. It will fall out from the 
eddy where it was at an earlier instant and will therefore rapidly decorrelate from the 
flow. In order to accommodate the back-flow necessary to satisfy continuity, the normal 
correlations must then contain negative loops (as those seen in figures]^ and d). Follow¬ 
ing Csanady (1963) we define the period of oscillation of the fluctuations as the ratio of 


the typical eddy diameter in the direction of gravity (i.e. the longitudinal integral scale 
Lq), and the particle terminal velocity Vt obtaining t = Lofvt — 16. This value is similar 
to the period of oscillations in the correlations in figure 

As shown in figure |^) , the quiescent environment presents a peculiar behavior of the 
settling velocity correlation at f = 1%. In particular we observe oscillations around zero 
of long period, T = 160 — 180. From the analysis of particle snapshots at different times 
(not shown), we observe that these seem correlated to the formation of regions of different 
density of the particle concentration. Hence a particle crossing regions with different local 
particle concentration may experience a varying settling velocity. 

To further understand the particle dynamics, we display in figure]^) and b) the single 
particle dispersion, i.e. the mean square displacement, for both quiescent and turbulent 
cases at ^ = 0.5% and 1.0%. The mean displacement, {Vp^r)t, is subtracted from the 
instantaneous displacement in the settling direction Az{t) to highlight the fluctuations 
with respect to the mean motion. For all cases we found initially a quadratic scaling in 
time {{Axf) ^ t^) typical of correlated motions while the linear diffusive behavior takes 
over at longer times ((Axf) ^ 2Det with He the diffusion coefficient). 

The turbulent cases show a similar behavior for both volume fractions. The crossover 
time when the initial quadratic scaling is lost and the linear one takes place is about 
t ~ 10 and t ~ 50 for the normal and tangential component, respectively. This difference 
is consistent with the correlation timescales previously discussed. The dispersion rates are 
similar in all directions in a turbulent environment. The quiescent cases present different 
features. First of all, dispersion is much more effective in the settling direction than in 
the normal one. The dispersion rate is smaller than the turbulent cases in the horizontal 
directions, while, surprisingly, the mean square displacement in the settling direction is 
similar to that of the turbulent cases being even higher at (f = 0.5%, something we relate 
to the drafting-kissing-tumbling behaviour discussed above. The crossover time scale is 
similar to that of turbulent cases with the exception of the most dilute case which does 
not reach a fully diffusive behavior at t ~ 500. This long correlation time makes the 
mean square displacement of this case higher than the corresponding turbulent case at 
long times. 
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Figure 7. Mean square particle displacement in the directions parallel and perpendicular to 
gravity for both quiescent and turbulent cases for a) = 0.5% and h) 4> = !%■ The mean particle 
displacement {Vp^T)t in the settling direction is subtracted from the instantaneous displacement 
when computing the statistics. In the inset, we show (Aa;f)/t as a function of time as well as 
the diffusion coefficients for the turbulent cases. 


In the insets of figure 13^) and b) we show {Axf)/t as a function of time. In all cases 
except the quiescent one at (/> = 0.5%, the diffusive regime is reached and it is possible to 
calculate the diffusion coefficients De = (Axf) /{2t). For the turbulent case with </> = 0.5% 
we obtain Dg = 0.52 and 0.28 in the directions parallel and perpendicular to gravity 
while for (p = 1% we obtain = 0.57 and 0.32. In the quiescent cases the diffusion 
coefficients in the horizontal directions are approximately 0.03, whereas the coefficient 
in the gravity direction a,t p = 1% is about 0.40. Csanady (19631 proposed a theoretical 
estimate of the diffusion coefficients for pointwise particles. Using these estimates, we 
obtain approximately Dg = 1.4 and 0.7 in the directions parallel and perpendicular to 
gravity. These are about 2.5 times larger than those found here for finite-size particles. 
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Quiescent f = 0.5% Turbulent f = 0.5% Quiescent f = 1% Turbulent <f = 1 
+0.18 +0.28 +0.25 +0.29 

avf[„ +0.04 +0.27 +0.06 +0.27 

Table 3. Fluctuation rms of the fluid velocities parallel avf ^ and perpendicular avf „ to 
gravity. The turbulent fluid velocity undisturbed rms is ~ 0.3 


3.2. Fluid statistics 

Table reports the fluctuation intensities of the fluid velocities for all cases considered. 
These are calculated by excluding the volume occupied by the spheres at each time step 
and averaging over the number of samples associated with the fluid phase volume. As 
expected, the fluid velocity fluctuations are smaller in the quiescent cases than in the 
turbulent regime. In the quiescent case the rms of the velocity fluctuations is about 50% 
larger In the settling direction than in the normal direction because of the long range 
disturbance induced by the particles wakes. The increase of the volume fraction enhances 
the fluctuations in both directions. Fluctuations are always larger in the turbulent case, 
with the most significant differences compared to the quiescent cases in the normal direc¬ 
tion, where the presence of the buoyant solid phase brakes the Isotropy of the turbulent 
velocity fluctuations. 

Hence, the solid phase clearly affects the turbulent flow field. Although the present 
study focuses on the settling dynamics, it is interesting to briefly discuss how turbulence is 
modulated. Modulation of isotropic turbulence by neutrally buoyant particles is examined 


Lucci et al. (2010); however the results change due to buoyancy as investigated here. 


Typical turbulent quantities are reported in table where they are compared with the 
unladen case at 0 = 0. The energy dissipation e increases with (j) becoming almost double 
at (j) = 1%. This behavior is expected since the buoyant particles inject energy in the 
system that is transformed into kinetic energy of the fluid phase that has to be dissipated. 
The higher energy flux, i.e. dissipation, is reflected in a reduction of the Kolmogorov 
length 77 . The particles reduce the velocity fluctuations, decreasing the turbulent kinetic 
energy level. The combined effect on k and e result in a decrease of the Taylor microscale 
A and of Rey, likewise the integral length Lq and Rclo also decrease. The reduction 
of the large and small turbulence scales is associated to the additional energy injection 
from the settling particles. Energy injection occurs at the size of the particles, which is 
below the unperturbed integral scale Lq explaining the lowering of the effective integral 
Lq and of Taylor A length-scales. This additional energy is transferred to the bulk flow 
in the particle wake. Associated to this energy input there is a new mechanism for 
dissipation that is the interaction of the flow with the no-slip surface of the particles. 
The mean energy dissipation held in the particle reference frame for the turbulent case 
with (j) = 0.5% is therefore shown in hgure[^ After a statistically steady state is reached, 
the norm of the symmetric part of the velocity gradient tensor Ey and the dissipation 
e = 2vEijEij are calculated at each time step on a cubic mesh centered around each 
particle; the dissipation is calculated on the grid points outside the particle volume. The 
data presented have been averaged over all particles and time to get the mean dissipation 
held displayed in the hgure. The maximum (e) is found around the particle surface with 
maximum values in the front; the mean dissipation drops down to the values found in 
the rest of the domain on the particle rear. The overall energy dissipation is therefore 
made up of two parts: the hrst associated to the dissipative eddies far from the particle 
surfaces and the second associated to the mean and huctuating how held near the particle 
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n/{2a) 

k 

A/(2a) 

Re\ 

e 

Te 

Rclo 

0.000 

0.084 

0.13 

1.56 

90 

0.0026 

47.86 

1205 

0.005 

0.077 

0.10 

1.19 

62 

0.0037 

27.97 

570 

0.010 

0.069 

0.11 

1.01 

54 

0.0055 

19.88 

435 


Table 4. Turbulent flow parameters in particle units for (j) = 0, <j) = 0.5% and (p = 1%. 



- 2-1012 

z/(2a) 


Figure 8. Contours of the mean fluid dissipation, (e) {vt /{2a)), averaged in the particle frame 

of reference. 

surface. To conclude, the settling strongly alters the typical turbulence features via an 
anisotropic energy injection and dissipation, thus breaking the isotropy of the unladen 
turbulent flow. The energy is injected by the fluctuations in the particle wake whereas 
stronger energy dissipation occurs in the front of each particle. As a consequence, the 
fluid velocity fluctuations change in the directions parallel and perpendicular to gravity 
as shown in table ID 


3.3. Relative velocity 

An important quantity to understand and model the settling dynamics is the particle to 
fluid relative motion. Although it is still unclear how to properly calculate the slip velocity 
between the two phases, we consider spherical shells around each particle, centered on 


the particles centroids, inspired by the works of Bellani & Variano (2012) and Cisse et al. 


(2013). We calculate the mean difference between the particle and fluid velocities in each 


shell as 


rel)x.,t,NP — ( Up 


1 


n{A) 


UfdV 


/ t,NP 

where r2(A) is the volume of a shell of inner radius A. A parametric study on the slip 
velocity is performed by changing the inner radii of these spherical shells from A = 0.75 
particle diameters to A = 5.0 particle diameters, while keeping the shell thickness 6 
constant and equal to 0.063 in units of 2a. In figure]^ we report the component of TJrei 
parallel to gravity as a function of the shells inner radii A/(2a). As the shell inner radii 
increase, \Urei,T\ tends exponentially to an asymptotic value which corresponds to the 
mean particle velocity in the same direction, {Vp^r)- This is expected since the correlation 
between the fluid and particle velocity goes to zero at large distances. The quiescent cases 
still show a 0.5% difference between \Urei,T\ and (Ijj.t) at A/(2a) = 5. This difference is 
again attributed to the long coherent wakes of the particles. 














Sedimentation of finite-size spheres 


19 



Figure 9. The component of Urei aligned with gravity as a function of the shell inner radii 
A/(2a) for the four cases studied. The dashed lines represent instead the mean particles velocities 
in the direction of gravity, 


A/(2a) 





0.75 

+0.81 

0.06 

+8.040 

123.29 

1.00 

+0.88 

0.09 

+5.655 

68.26 

1.25 

+0.92 

0.10 

+4.547 

47.66 

1.50 

+0.93 

0.11 

+3.912 

37.69 

1.75 

+0.94 

0.12 

+3.482 

31.83 

2.50 

+0.95 

0.13 

+2.919 

24.34 

4.00 

+0.95 

0.14 

+2.18 

16.25 

5.00 

+0.96 

0.14 

+2.22 

17.09 


Table 5. Moments of the pd/(?7rei,r) for f = 0.5% in the quiescent case. The thickness of the 
shell used to compute the slip velocity is S/{2a) — 0.063. 


A/(2a) 

{Urel ,t'} 




0.75 

+0.76 

0.07 

+0.072 

3.90 

1.00 

+0.83 

0.09 

+0.010 

4.28 

1.25 

+0.85 

0.10 

+0.087 

4.56 

1.50 

+0.86 

0.11 

+0.128 

4.54 

1.75 

+0.87 

0.11 

+0.117 

4.32 

2.50 

+0.87 

0.13 

+0.054 

3.91 

4.00 

+0.88 

0.16 

+0.046 

3.32 

5.00 

+0.88 

0.18 

+0.047 

3.10 


Table 6. Moments of the pdf{Urei,T) for f = 0.5% in a turbulent environment. The thickness 
of the shell used to compute the slip velocity is S/{2a) — 0.063. 


The probability density function of these relative velocities, Urel,T^ and their first four 
central moments are computed and reported in tables and as a function of A/(2a) 
for the quiescent fluid and turbulent cases at (^ = 0.5%. In the turbulent case, the 
moments approach those of a Gaussian distribution with vanishing skewness and flatness 
close to 3, especially at large A. In the quiescent case, the third and fourth moments 
display higher values that decrease as A is increased, tending to the values of the particle 
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Figure 10. Comparison of the probability density functions of tangential relative velocity and 
corresponding particle Reynolds number Rcp. Colors and symbols are the same as previous 
figures. In Panel a) the cases for A/(2a) = 1.5 are reported while in panel b) we show the 
results for A/(2a) = 2.5. 


velocities. The pdfs pertaining the four cases considered, calculated in spherical shells 
with an inner radius of 1.5 and 2.5 particle diameters, are compared in figure |10[ A 
second axis, reporting the the particle Reynolds number Rcp based on Urei,T, is also 
displayed in each hgure. In the former case, A/(2a) = 1.5, the shell radius is of the 
order of the Taylor scale to highlight the particle dynamics, while the relative velocity is 
approaching the asymptotic values for the larger shell. The pdfs of the relative velocity 
appear narrower than those of the particle absolute velocity, indicating that the particles 
tend to be transported by the large-scale motions, filtering the smallest scales. 

The distributions pertaining the simulations in a turbulent environment are nearly 
Gaussian with modal values well below one. The quiescent cases show skewed distribu¬ 
tions with long tails at high velocities, as observed for the particle velocities in figure]^). 
The particles settle on average with a velocity close to that of a single particle, with occa¬ 
sional events of higher velocity due to the drafting-kissing-tumbling dynamics. The lower 
the volume fraction the more intermittent is the dynamics. 

Knowing the tangential fluid velocity {Vf^r)r{p,t), averaged in each shell and at each 
time step, and the corresponding tangential particle velocities (p, t), it is then possible 
to find their joint probability distribution (for sake of simplicity we write 

{Vf,r)r as These are evaluated in shells at A/(2a) = 1.75 for each case studied and 
reported in figure 


11 


Vp,r), 


In each case, the integral of P{Vf^ 

/ OO 

yp,T)yf,T dyf^T, 

-OO 


(3.3) 


is also reported (continuous lines). This represents the most probable particle velocity 
Vp^T given a certain fluid velocity V/.r, or, equivalently, the most probable fluid velocity 
surrounding a particle settling with velocity In the turbulent cases these integrals 
are well approximated by straight lines (displayed with dashed lines in the hgure) 




rp,r=Ciyf,r+C2. (3.4) 

In both cases Ci is approximately 1 while C 2 is about 0.86 for (f = 0.5% and 0.84 for 
(j) = 1%. These values are in agreement with the values found for the average relative 
velocities of shells at A/(2o) = 1.75. In a quiescent how, conversely, the integral in 
eq. (3.31 gives a curved line and no best-ht is therefore reported. In these cases, the 
joint probability distribution is broader, particularly in the region of higher particle 
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Figure 11. Joint probability distributions of Vf,T and V^,r evaluated in spherical shells located 
at 1.75 particle diameters from each particle. We report in panels a) and c) the quiescent cases 
for f = 0.5% and 1%, while the respective turbulent cases are reported in panels b) and d). The 
continuous lines represent the integrals of the joint probability distributions. In the turbulent 
cases the dashed lines represent the best-£t of these integrals. 


velocities, Vp^r- This is again due to the intense particle interactions and the drafting- 
kissing-tumbling behaviour described in figure]^ which confirms the high flatness of the 
probability density functions of the relative particle velocities. 

Further insight can be obtained by plotting isocontours of the average particle relative 
velocities and their fluctuations, in both quiescent and turbulent flows. To this end, we 
follow the approach by Garcia-Villalba et al. (2012). We place an uniform and structured 
rectangular mesh around each particle, with origin at the particle center. By means of 
trilinear interpolations we And the fluid and relative velocities on this local mesh and 
average over time and the number of particles to obtain the mean relative velocity held 
and its fluctuations. 

The vertical velocity of the single sphere settling in a quiescent fluid is reported in 
figure [T^ The relative normal and tangential velocities and their fluctuations are instead 
displayed in figure 13 for suspensions with (f = 0.5% in both quiescent and turbulent 
environment. In the quiescent fluid simulations, a long wake forms behind the represen¬ 
tative particle and, as seen from the single point particle velocity correlations, it takes a 
long time for this velocity fluctuations to decorrelate. In the turbulent case instead, the 
wakes are disrupted by the background fluctuations. 

Interesting observations can be drawn from the relative velocity fluctuation fields. 
Comparing figures [Ts] :) and 131) we note that intense vortex shedding occurs around the 
particles in the turbulent case, with important fluctuations of Urei,T- From figures 13 3) 
















22 


W. Fornari, F. Picano and L. Brandt 



0.2 


-6 -4-2 0 2 

z/(2a) 

Figure 12. Contour plot of the velocity component in the direction of gravity for the single 

sphere settling in a quiescent fluid. 


and 13') we also see that the relative velocity fluctuations are drastically increased in the 


horizontal directions in a turbulent environment. Noteworthy, vortex shedding occurs at 
particle Reynolds numbers below the critical value above which this is usually observed 
(Bouchet et al. 20061. Vortex shedding is unsteady in nature and unsteady effects may 
therefore play an important role in the increase of the overall drag, as further discussed 
below. Lower fluctuation intensities are found on the front part of the particles, where 
the energy dissipation is highest, and the immediate wake in the recirculating region 
where the instability is found to develop. 


3.4. Drag analysis 


As in Maxey & Riley (1983), we write the balance of the forces acting on a single sphere 
settling through a turbulent flow. The equation of motion for a spherical particle reads 

4 , dV„ 


-TTa^Pp- 


dt 


= ^™^(Pp-P/)g+ 

JdV„ 


ndS 


(3.5) 


where the integral is over the surface of the sphere dVp, n is the outward normal and 
T = —pi + 2fiE is the fluid stress. As commonly done in aerodynamics, we replace the 
last term of equation (3.5) with a term depending on the relative velocity Vrei and a 
drag coefficient Cd 




—TTa 

3 


Pp- 


dt 


= ^7ra^(pp - Pf)g - ^Pfira^lUrell’^relCo 


(3.6) 


with 7ra^ the reference area. Generally, the drag coefficient Cd is a function of a Reynolds 
number and a Strouhal number which accounts for unsteady effects. In the present case, 
we consider it to be a function of the Reynolds number based on the relative velocity 
Rcp = 2a|Ur.ez|/u (in a turbulent held it is proper to define Rcp in terms of the relative 
velocity between particle and fluid) and of the Strouhal number defined as 


St = 


^( 2 «) 

lU.ezP ■ 


(3.7) 


The drag on the particle depends on both nonlinear and unsteady effects (such as the 
Basset history force and the added mass contribution) through these two non-dimensional 
numbers. 

Commonly, the unsteady contribution is neglected and Cr? is assumed to depend only 
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Figure 13. Fields of {[/rei.r), Urei,Ti 71 for the quiescent (left column) and turbulent cases 

(right column). 


on the Reynolds number. Since we want to investigate both nonlinear and unsteady ef¬ 
fects, we decide to express the drag coefficient as CD{Rep, St) = Cuo {R^p) [1 + f^iR^p, S't)], 
yielding 


dt 


= ^7ra^(Pp - p/)g - ^Pfna'^\'Urei\'VreiCDo{Rep) [1 + i^iRep, St)], (3.8) 


where ■i/i = 0 V Rcp if 5”^ = 0 (steady motion). We can therefore identify a quasi-steady 
term and the extra term which accounts for unsteady phenomena. 

By ensemble averaging equation (3.8) over time and the number of particles, and 
assuming the settling process to be at statistically steady state, we can find the most 
important contributions to the overall drag. The steady-state average equation reads 


1 


0 = -Tra-^{pp - pf)g - -Pfira^ {\lJrei\'UreiCDoiRep) [1 + ifiRep, St)]). 


(3.9) 


Denoting the entire time dependent term simply as 'F(t) and rearranging, we obtain the 
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following balance 
4 


-na 

3 


Pp y 
Pf 


g = -7ra'^{\lJrel\'^relCDoiRep)) + 4'(t). 


(3.10) 


The term on the left-hand side is known whereas the time-dependent term 4'(t) is of 
difficult evaluation. The nonlinear steady term can be calculated using the approach de¬ 
scribed in section i2. At each time step we calculate the relative velocity in the spherical 
shells surrounding each particle. From these we compute the particle Reynolds number 
Rcp = 2a\Urei\/v and using equation (2.8) (where we replace the terminal Reynolds 


number with the new particle Reynolds number) we obtain the drag coefficient. The first 
term on the right-hand side can therefore be evaluated by averaging over the number 
of particles and time steps considered. Finally we divide everything by the buoyancy 
acceleration to estimate the relative importance of each term 


100% = S{Rep) -b 4'(t) 


(3.11) 


where S{Rep) represents the nonlinear steady term while the unsteady term has been 
denoted again as 4'(t) for simplicity. 

This approach is applied to the results of the simulations of a single sphere and to the 
quiescent and turbulent cases with (j) = 0.5% and 1%. The inner radius of the sampling 
shells is chosen to be 5 particle radii and the results obtained are reported in table 
The single sphere simulation provides an estimate of the error of the method. Since our 
terminal Reynolds number is smaller than the critical Reynolds number above which 
unsteady effects become important and our velocity field is indeed steady, the term 4' (t) 
should be negligible. The critical Reynolds number Rcc has been found to be approxi¬ 


mately 274 by Bouchet et al. (2006) while we recall that our terminal Reynolds number 


for the isolated particle is 200. The nonlinear steady term provides in this case, however, 
an overestimated value of the drag, with a percentage error of about -1-3% with respect 
to the buoyancy term. The possible causes of this error are the long wake and the fact 


that equation (2.8) is empirical. The data from the single-particle simulations are used 
to correct the results from the other runs, i.e. the data are normalized with the total 
drag obtained in this case. In the quiescent case at ^ = 0.5%, unsteady effects are neg¬ 
ligible (about 0.5% of the total drag), while they increase to about 6% when increasing 
the particle evolve fraction to 1%. In a turbulent flow, importantly, we notice that the 
contribution of 4'(t) adds up to approximately 10% of the total at <() = 0.5% and to about 
12 % of the total at (p = 1%. 

Note that one can write the steady drag as mean and fluctuating component 

SiRep) = {\Vrel\VrelCDoiRep)} = {Urelf C Doi{Rep)) + S\{Rep)). 

The fluctuations S'{{Rep)) would be responsible of the reduction of the settling velocity 
if this were to be attributed to nonlinear drag effects only (see also Wang & Maxey 
1993). We verified that for our results, the total and mean component differ by about 


2%, {\TJrei\^reiCDo{Rep)) « {Urei)"^CDo{{Re-p))■ This leads us to the conclusion that 
the main contribution to the overall drag is due to the steady term but the reduction 
of the mean settling velocity in a turbulent environment is almost entirely due to the 
various unsteady effects. These can be related to unsteady vortex shedding, see figure [T^ 
as in the experiments of a single sphere by Mordant & Pinton (2000) These observations 


are also in agreements with the results in Homann et al. (2013). These authors observe 


that the enhancement of the drag of a sphere towed in a turbulentt environment can be 
explained by the modification of the mean velocity and pressure profile, and thus of the 
boundary layer around the sphere, by the turbulent fluctuations. 
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case S{Rep) 

Quiescent f = 0.5% 99.5% 
Turbulent f = 0.5% 90.4% 
Quiescent f = 1.0% 94.1% 
Turbulent f = 1.0% 87.8% 


4 /( 1 ) 

0.5% 

9.6% 

5.9% 

12 . 2 % 


Table 7. Percentage contribution of the non-linear and unsteady terms for the quiescent and 
turbulent case with f = 0.5%. The data are normalized by the mean drag from the simulation 
of a single sphere in quiescent fluid. 


4. Final remarks 

We report numerical simulations of a suspension of rigid spherical slightly-heavy par¬ 
ticles in a quiescent and turbulent environment using a direct-forcing immersed bound¬ 
ary method to capture the fluid-structure interactions. Two dilute volume fractions, 
(j) = 0.5% and 1%, are investigated in quiescent fluid and homogeneous isotropic turbu¬ 
lence at Re\ = 90. The particle diameter is of the order of the Taylor length scale and 
about 12 times the dissipative Kolmogorov scale. The ratio between the sedimentation 
velocity and the turbulent fluctuations is about 3.4, so that the strongest fluid-particle 
interactions occur at approximately the Taylor scale. 

The choice of the parameters is inspired by the reduction in sedimentation velocity 


observed experimentally in a turbulent flow by Byron (20151 and in the group of Prof. 


Variano at UC Berkeley. In the experiment, the isotropic homogeneous turbulence is 
generated in a tank of dimensions of several integral lengthscales by means of two facing 
randomly-actuated synthetic jet arrays (driven stochastically). The Taylor microscale 
Reynolds number of the experiment is Rex = 260. Particle Image Velocimetry using 
Refractive-Index-Matched Hydrogel particles is used to measure the fluid velocity and 
the linear and angular velocities of finite-size particles of diameter of about 1.4 Taylor 
microscales and density ratios Pp/p/ = 1.02, 1.006 and 1.003. The ratio between the 
terminal quiescent settling velocity Vt and the turbulence fluctuating velocity u' is about 


1, higher than in our simulations where this ratio is 3.3. Byron (20151 observes reductions 


of the slip velocity between 40% and 60 % when varying the shape and density of the 


particles. As suggested by Byron (2015), the larger reduction in settling velocity observed 


in the experiments is most likely explained by the larger turbulence intensity. 

The new findings reported here can be summarized as follows: i) the reduction of set¬ 
tling velocity in a quiescent flow due to the hindering effect is reduced at finite inertia 
by pair-interactions, e.g. drafting-kissing-tumbling. ii) Owing to these particle-particle 
interactions, sedimentation in quiescent environment presents therefore significant inter- 
mittency. iii) The particle settling velocity is further reduced in a turbulent environment 
due to unsteady drag effects, iv) Vortex shedding and wake disruption is served also in 
subcritical conditions in an already turbulent flow. 

In a quiescent environment, the mean settling velocity slightly decreases from the 
reference value pertaining few isolated particle when the volume fraction </> = 0.5% and 
(j) = 1%. This limited reduction of the settling velocity with the volume fraction is in 
agreement with previous experimental findings in inertia-less and inertial flows. The 
Archimedes number of our particle is 21000, in the steady vertical regime before the 
occurrence of a first bifurcation to an asymmetric wake. In this regime, Uhlmann & 


Doychev (2014) observe no significant particle clustering, which is is confirmed by the 


present data. 

The skewness and flatness of the particle velocity reveal large positive values in a 
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quiescent fluid, and accordingly the velocity probability distribution functions display 
evident positive tails. This indicates a highly intermittent behavior. In particular, it is 
most likely to see particles sedimenting at velocity significantly higher than the mean: 
this is caused by the close interactions between particle pairs (more seldom triplets). 
Particles approaching each other draft-kiss-tumble while falling faster than the average. 

In a turbulent flow, the mean sedimentation velocity further reduces, to 0.88 and 0.86 
at (j) = 0.5% and (j) = 1%. The variance of both the linear and angular velocity increases in 
a turbulent environment and the single-particle time correlations decay faster due to the 
turbulence mixing. The velocity probability distribution function are almost symmetric 
and tend towards a Gaussian of corresponding variance. The particle lateral dispersion is, 
as expected, higher in a turbulent flow, whereas the vertical one is, surprisingly, of com¬ 
parable magnitude in the two regimes; this can be explained by the highly intermittent 
behavior observed in the quiescent fluid. 

We compute the averaged relative velocity in the particle reference frame and the 
fluctuations around the mean. We show that the wake behind each particle is on average 
significantly reduced in the turbulent flow and large-amplitude unsteady motions appear 
on the side of the sphere in the regions of minimum pressure where vortex shedding is 
typically observed. The effect of a turbulent flow on the damping of the wake behind 
a rigid sphere has been discussed for example by Bagchi & Balachandar (2003), while 
the case of a spherical bubble has been investigated by Merle et al. (2005). Using the 
slip velocity between the particle and the fluid surrounding it, we estimate the nonlinear 
drag on each particle from empirical formulas and quantify the relevance of non-stationary 
effects on the particle sedimentation. We show that these become relevant in the turbulent 
regime, amount to about 10-12% of the total drag, and are responsible for the reduction 
of settling velocity with the respect to the quiescent flow. This can be compared with 
the simulations in Good et al. (2014) who attribute the reduction of the sedimentation 
velocity of small (2a < 77 ) heavy {pp/pf « 900) spherical particles in turbulence to the 
nonlinear drag. Here, we show that non-stationary effects become relevant for larger 
particles at lower density ratios. 

The present investigation can be extended in a number of interesting directions. Pre¬ 
liminary simulations reveal that variations of the density ratio at constant Archimedes 
number do not significantly modify the results presented here. It would be therefore in¬ 
teresting to investigate the effect of the Galileo number on the particle dynamics and of 
the ratio between turbulent fluctuations and sedimentation velocity. 
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